
using DataFrames, Plots, Statistics, JSON, Distributed
plotlyjs();


const numcores = 31

if nprocs() < numcores
    addprocs(numcores - nprocs())
end

@show nprocs();

nprocs() = 31


@everywhere include(joinpath(dirname(pwd()),"jgrc/TuLiPa/src/TuLiPa.jl"));


@everywhere include(joinpath(dirname(pwd()),"jgrc/JulES/src/JulES.jl"));


datayear = 2025
scenarioyearstart = 1981
scenarioyearstop = 2011
totalscen = 30 # scenarios to consider uncertainty for

# Standard time for market clearing - perfect information so simple time type
tnormal = FixedDataTwoTime(getisoyearstart(datayear),getisoyearstart(scenarioyearstart))

# Phasein settings
phaseinoffsetdays = 2 # also simulation step length
phaseinoffset = Millisecond(Day(phaseinoffsetdays)) # phase in straight away from second stage scenarios
phaseindelta = Millisecond(Week(5)) # Phase in the second stage scenario over 5 weeks
phaseinsteps = 5 # Phase in second stage scenario in 5 steps

# Make scenario times for all uncertainty scenarios. List of tuples with tnormal, tphasein and scenarionumber
totalscentimes = []
for scen in 1:totalscen
    scentnormal = FixedDataTwoTime(getisoyearstart(datayear),getisoyearstart(scenarioyearstart + scen - 1))
    scentphasein = PhaseinTwoTime(getisoyearstart(datayear),getisoyearstart(scenarioyearstart), getisoyearstart(scenarioyearstart + scen - 1), phaseinoffset, phaseindelta, phaseinsteps);

    push!(totalscentimes, (scentnormal, scentphasein, scen))
end


sti_themadata = "data_fra_thema"

themastructure = json_to_elements(sti_themadata, "dataset_thema.json")
themaseries = json_to_elements(sti_themadata, "tidsserier_thema.json")

elements = vcat(themaseries,themastructure)
addscenariotimeperiod!(elements, scenarioyearstart, scenarioyearstop);


# Set horizons for price prognosis models
# Long
longhorizonduration = Millisecond(Week(5*52))
longhydroperiodduration = Millisecond(Week(6))

longrhsdata = StaticRHSAHData("Power", datayear, scenarioyearstart, scenarioyearstop)
longmethod = KMeansAHMethod()
longclusters = 4
longunitduration = Millisecond(Hour(6))

# Medium
medhorizonduration = Millisecond(Week(54))
medhydroperiodduration = Millisecond(Day(7)); @assert medhorizonduration.value % longhydroperiodduration.value == 0

medrhsdata = DynamicRHSAHData("Power")
medmethod = KMeansAHMethod()
medclusters = 4
medunitduration = Millisecond(Hour(4))

# Short
shorthorizonduration = Millisecond(Week(1))
shorthydroperiodduration = Millisecond(Day(1)); @assert medhorizonduration.value % shorthorizonduration.value == 0
shortpowerparts = 12

# Preallocate storage for problems and results on different cores. Use package DistributedArrays
# Distribute scenarios
allscenarios = distribute(totalscentimes)

# Problems are built, updated, solved, and stored on a specific core. Moving a problem between cores is expensive, so we want it to only exist on one core. 
longprobs = distribute([HiGHS_Prob() for i in 1:length(allscenarios)], allscenarios)
medprobs = distribute([HiGHS_Prob() for i in 1:length(allscenarios)], allscenarios)
shortprobs = distribute([HiGHS_Prob() for i in 1:length(allscenarios)], allscenarios)

# Results are moved between cores. These are much smaller than longprobs/medprobs/shortprobs and are inexpensive to move between cores.
medprices = distribute([Dict() for i in 1:length(allscenarios)], allscenarios)
shortprices = distribute([Dict() for i in 1:length(allscenarios)], allscenarios)
medendvaluesobjs = distribute([EndValues() for i in 1:length(allscenarios)], allscenarios)
nonstoragestates = distribute([Dict{StateVariableInfo, Float64}() for i in 1:length(allscenarios)], allscenarios)

# Organise inputs and outputs
probs = (longprobs, medprobs, shortprobs)
longinput = (longhorizonduration, longhydroperiodduration, longrhsdata, longmethod, longclusters, longunitduration, "long")
medinput = (medhorizonduration, medhydroperiodduration, medrhsdata, medmethod, medclusters, medunitduration, "long")
shortinput = (shorthorizonduration, shorthydroperiodduration, shortpowerparts, "short")
allinput = (numcores, elements, allscenarios, phaseinoffset)
output = (medprices, shortprices, medendvaluesobjs, nonstoragestates)
  
# Initialize price prognosis models and run for first time step. Run scenarios in parallell
@time pl_prognosis_init!(probs, allinput, longinput, medinput, shortinput, output)

 86.040697 seconds (5.53 M allocations: 252.401 MiB, 1.22% compilation time)

Task (done) @0x00000000744e3650


scenario = 1
index = collect(medprices[scenario]["steprange"])
prices = medprices[scenario]["matrix"]
labels = [name for name in medprices[scenario]["names"]]
p = plot(index,prices,label=reshape(labels, 1, length(labels)),legend=:outertopright)

for scenario in 2:5 # length(allscenarios)
    prices = medprices[scenario]["matrix"]
    labels = [name for name in medprices[scenario]["names"]]
    plot!(p,index,prices,label=reshape(labels, 1, length(labels)),legend=:outertopright)
end
display(p)


scenario = 1
index = collect(shortprices[scenario]["steprange"])
prices = shortprices[scenario]["matrix"]
labels = [name for name in shortprices[scenario]["names"]]
p = plot(index,prices,label=reshape(labels, 1, length(labels)), ylabel="€/MWh",legend=:outertopright)

for scenario in 2:5 # length(allscenarios)
    prices = shortprices[scenario]["matrix"]
    labels = [name for name in shortprices[scenario]["names"]]
    plot!(p, index,prices ,label=reshape(labels, 1, length(labels)),legend=:outertopright)
end
display(p)


# Read dataset for stochastic and market clearing models
sti_themadata = "data_fra_thema"
tst = JSON.parsefile(joinpath(sti_themadata, "dataset_thema_excl_hydro_nose.json")) # European dataset without aggregated hydro in Norway and Sweden
themastructure = getelements(tst) 
tse = JSON.parsefile(joinpath(sti_themadata, "tidsserier_thema.json"))
themaseries = getelements(tse, sti_themadata)

dse = JSON.parsefile("data_fra_dynmodell/tidsserier_detd.json")
detdseries = getelements(dse);
dst = JSON.parsefile("data_fra_dynmodell/dataset_detd.json")
detdstructure = getelements(dst);

detailedelements = vcat(themaseries,themastructure,detdseries,detdstructure)
addscenariotimeperiod!(detailedelements, scenarioyearstart, scenarioyearstop);


# Mapping between aggregated and detailed storages
detailedrescopl = JSON.parsefile("data_fra_dynmodell/magasin_elspot1.json")

# Global energy equivalent detailed reservoirs
enekvglobaldict = Dict()
for element in detailedelements
    if element.typename == GLOBALENEQKEY
        enekvglobaldict[split(element.instancename,"GlobalEneq_")[2]] = element.value["Value"]
    end
end

# Detailed dataset has reservoirs for SE4, aggregated does not, TODO: Improve aggregation/mapping
for k in keys(detailedrescopl)
    if detailedrescopl[k] == "SE4"
        detailedrescopl[k] = "SE3"
    end
end

# Get dictionary with each detailed reservoir and their water value for each scenario
# TODO: Detailed run-of-river reservoirs get water value from aggregated reservoir hydro
function getendvaluesdicts(endvaluesobjs, detailedrescopl::Dict, enekvglobaldict::Dict)
    endvaluesdicts = Dict[]
    for endvaluesobj in endvaluesobjs
        instance = [getinstancename(getid(obj)) for obj in endvaluesobj.objects]
        endvalues = endvaluesobj.values
        endvaluesdict = Dict(instance .=> endvalues)

        for (k,v) in detailedrescopl
            endvaluesdict["Reservoir_" * k] = endvaluesdict["Reservoir_" * v * "_hydro_reservoir"] * enekvglobaldict[k]
        end
        push!(endvaluesdicts, endvaluesdict)
    end
    
    return endvaluesdicts
end
medendvaluesdicts = getendvaluesdicts(medendvaluesobjs, detailedrescopl, enekvglobaldict);


# Scenario reduction to this amount
numscen = 7

# Modelobjects that can be used to reduce scenarios
scenarioelements = copy(detailedelements)

# Horizons are needed to build modelobjects, but not used in scenario modelling
dummyperiods = 10
dummyperiodduration = Millisecond(Hour(24))
power_horizon = SequentialHorizon(dummyperiods, dummyperiodduration)
hydro_horizon = SequentialHorizon(dummyperiods, dummyperiodduration)

set_horizon!(scenarioelements, "Power", power_horizon)
set_horizon!(scenarioelements, "Battery", power_horizon)
set_horizon!(scenarioelements, "Hydro", hydro_horizon)

scenarioobjects = collect(values(getmodelobjects(scenarioelements)))

# Scenario modelling method
scendelta = MsTimeDelta(Day(364)) # scenario modelling based on the next year, even though the scenario problems can be longer
parts = 4 # divide scendelta into this many parts, calculate sum inflow for each part of the inflow series, then use clustering algorithm
scenmodmethod = InflowClusteringMethod(numscen, parts)
# scenmodmethod = NoScenarioModellingMethod(numscen, totalscentimes)
@time scenariomodelling!(scenmodmethod, scenarioobjects, numscen, totalscentimes, scendelta); # see JulES/scenariomodelling.jl

  5.455343 seconds (15.31 M allocations: 699.816 MiB, 3.38% gc time, 96.53% compilation time)


# Cut parameters
maxcuts = 13 # preallocate fixed number of cuts, no cut selection
lb = -1e5 # lower bound of the future value in the first iteration
reltol = 0.0001 # relative tolerance

# Parameters for stochastic subsystem problems (could also split totalduration into master- and subduration)
shortstartstorage = 50
medstartstorage = 65
smpdp = Millisecond(Hour(2)) # short/med - master/sub - period duration - power/hydro (commodity)
smpdh = Millisecond(Hour(2))
sspdp = Millisecond(Hour(2))
sspdh = Millisecond(Hour(2)) # both master and subproblems for PHS and batteries has 2 hour resolution
mmpdp = Millisecond(Hour(24))
mmpdh = Millisecond(Hour(24)) # daily resolution in hydro master problems
mspdp = Millisecond(Hour(168))
mspdh = Millisecond(Hour(168)) # weekly resolution in hydro subproblems
shorttotalduration = shorthorizonduration # total duration of master and subproblem
medtotalduration = medhorizonduration - Millisecond(Week(2)) # we reuse prices for two weeks, so have to be two weeks shorter than price prognosis problem

# Make sure time resolution of hydro and power are compatible (TODO: Could add function that makes them compatible)
@assert ceil(Int64, phaseinoffset/smpdp) == ceil(Int64, phaseinoffset/smpdh)
@assert ceil(Int64, (shorttotalduration-phaseinoffset)/sspdp) == ceil(Int64, (shorttotalduration-phaseinoffset)/sspdh)
@assert ceil(Int64, phaseinoffset/mmpdp) == ceil(Int64, phaseinoffset/mmpdh)
@assert ceil(Int64, (medtotalduration-phaseinoffset)/mspdp) == ceil(Int64, (medtotalduration-phaseinoffset)/mspdh)

# Convert DistributedArray of prices to local process
medpriceslocal = convert(Vector{Dict}, medprices)
shortpriceslocal = convert(Vector{Dict}, shortprices)

# Inputs
storageinfo = (shortstartstorage, medstartstorage, medendvaluesdicts)
shortterminputs = (detailedelements, shorttotalduration, smpdp, smpdh, sspdp, sspdh, scenmodmethod.scentimes, phaseinoffset, shortpriceslocal, true)
medterminputs = (detailedelements, medtotalduration, mmpdp, mmpdh, mspdp, mspdh, scenmodmethod.scentimes, phaseinoffset, medpriceslocal, false)

ustoragesystemobjects = []
ushorts = []
# Make modelobjects for short-term subsystems
@time makemastersubobjects!(shortterminputs, ustoragesystemobjects, ushorts)
# Make modelobjects for medium-term subsystems
@time makemastersubobjects!(medterminputs, ustoragesystemobjects, ushorts)

# Distribute subsystems with inputs and outputs on different cores
storagesystemobjects, shorts = distribute_subsystems(ustoragesystemobjects, ushorts) # somewhat smart distribution of subsystems to cores based on how many modelobjects in eac subsystem
masters = distribute([HiGHS_Prob() for i in 1:length(storagesystemobjects)], storagesystemobjects)
subs = distribute([[] for i in 1:length(storagesystemobjects)], storagesystemobjects)
states = distribute([Dict{StateVariableInfo, Float64}() for i in 1:length(storagesystemobjects)], storagesystemobjects)
cuts = distribute([SimpleSingleCuts() for i in 1:length(storagesystemobjects)], storagesystemobjects)
storagesystems = distribute([Dict() for i in 1:length(storagesystemobjects)], storagesystemobjects)

# Initialize subsystem problems and run for first time step. Run subsystems in parallell
@time pl_stochastic_init!(numcores, storagesystemobjects, shorts, masters, subs, states, cuts, storageinfo, lb, maxcuts, reltol, scenmodmethod.scentimes)

  5.879978 seconds (16.36 M allocations: 706.884 MiB, 2.90% gc time, 52.02% compilation time)
  3.243251 seconds (12.81 M allocations: 518.070 MiB, 4.68% gc time, 2.44% compilation time)
  6.915521 seconds (2.88 M allocations: 100.548 MiB, 4.53% compilation time)

Task (done) @0x00000000735f22f0


# Bring data to local core
masterslocal = convert(Vector{Prob}, masters)
cutslocal = convert(Vector{SimpleSingleCuts}, cuts)
nonstoragestateslocal = convert(Vector{Dict}, nonstoragestates)

# Initialize market clearing problem and run for first time step
cpdp = Millisecond(Hour(2)) # clearing period duration power/battery
cnpp = ceil(Int64, phaseinoffset/cpdp) # clearing numperiods power/battery
cpdh = Millisecond(Hour(6)) # clearing period duration hydro
# cpdh = Millisecond(Hour(2)) # clearing period duration hydro
cnph = ceil(Int64, phaseinoffset/cpdh) # clearing numperiods hydro
@time clearing, nonstoragestatesmean, varendperiod = clearing_init(detailedelements, tnormal, phaseinoffset, cpdp, cpdh, masterslocal, cutslocal, nonstoragestateslocal);

 12.962708 seconds (16.55 M allocations: 990.515 MiB, 5.99% gc time, 53.05% compilation time)


# Initialize start states for next time step, also mapping to TheMA storages and max capacity in TheMA-storages
startstates = startstates_init(clearing, detailedrescopl, enekvglobaldict, longprobs[1], tnormal);


# Initialize and collect prices and start states
price = Dict()
powerhorizonix = argmax(getnumperiods(h) for h in clearing.horizons)
getareaprices!(price, clearing, clearing.horizons[powerhorizonix], tnormal)
areanames = price["names"]
ix = collect(price["steprange"])
pricematrix = copy(price["matrix"])

statenames = collect(keys(startstates))
statematrix = collect(values(startstates));


clearingobjects = Dict(zip([getid(obj) for obj in clearing.objects],clearing.objects)) # collect results from all areas
# resultobjects = getpowerobjects(clearingobjects,["NO2"]); # only collect results for one area
resultobjects = clearing.objects # collect results for all areas
prices, rhstermvalues, production, consumption, hydrolevels, batterylevels, powerbalances, rhsterms, rhstermbalances, plants, plantbalances, plantarrows, demands, demandbalances, demandarrows, hydrostorages, batterystorages = init_results(clearing, clearingobjects, resultobjects, cnpp, cnph, cpdp, tnormal, false);


# Only do scenario modelling and calculate new cuts every 8 days (other reuse scenarios and cuts)
skipmed = Millisecond(Day(0))
skipmax = Millisecond(Day(6))

# How many time steps to run the simulation for
steps = 16;
# steps = (getisoyearstart(scenarioyearstop) - getisoyearstart(scenarioyearstart)).value/phaseinoffset.value;


step = 2 # already ran first step in initialization

totaltime = @elapsed while step <= steps # while step <= steps and count elapsed time
    # Increment simulation/main scenario and uncertainty scenarios
    tnormal += phaseinoffset
    display(tnormal)

    for i in 1:length(totalscentimes)
        (scentnormal, scentphasein, scenario) = totalscentimes[i]
        scentnormal += phaseinoffset
        scentphasein = PhaseinTwoTime(getdatatime(scentnormal), getscenariotime(scentnormal), getscenariotime(scentnormal), phaseinoffset, phaseindelta, phaseinsteps)
        totalscentimes[i] = (scentnormal, scentphasein, scenario)
    end

    # Increment skipmed - should we reuse watervalues this time step?
    skipmed += Millisecond(phaseinoffset)
    if skipmed > skipmax
        skipmed = Millisecond(0)
    end

    # Deterministic long/mid/short - calculate scenarioprices for all 30 scenarios
    allscenarios = distribute(totalscentimes, allscenarios) # TODO: Find better solution
    @time pl_prognosis!(numcores, longprobs, medprobs, shortprobs, medprices, shortprices, nonstoragestates, startstates, allscenarios, skipmed)

    # Stochastic sub systems - calculate storage value    
    if skipmed.value == 0
        # Choose new scenarios
        scenariomodelling!(scenmodmethod, scenarioobjects, numscen, totalscentimes, scendelta)
        
        medpriceslocal = convert(Vector{Dict}, medprices)
        medendvaluesdicts = getendvaluesdicts(medendvaluesobjs, detailedrescopl, enekvglobaldict)
    else
        # Increment existing scenarios
        for i in 1:length(scenmodmethod.scentimes)
            (scentnormal, scentphasein, scenario) = scenmodmethod.scentimes[i]
            scentnormal += phaseinoffset
            scentphasein = PhaseinTwoTime(getdatatime(scentnormal), getscenariotime(scentnormal), getscenariotime(scentnormal), phaseinoffset, phaseindelta, phaseinsteps)
            scenmodmethod.scentimes[i] = (scentnormal, scentphasein, scenario)
        end
    end
    shortpriceslocal = convert(Vector{Dict}, shortprices)

    @time pl_stochastic!(numcores, masters, subs, states, cuts, startstates, medpriceslocal, shortpriceslocal, medendvaluesdicts, shorts, reltol, scenmodmethod.scentimes, skipmed)

    # Market clearing
    masterslocal = convert(Vector{Prob}, masters)
    cutslocal = convert(Vector{SimpleSingleCuts}, cuts)
    nonstoragestateslocal = convert(Vector{Dict}, nonstoragestates)

    @time clearing!(clearing, tnormal, startstates, masterslocal, cutslocal, nonstoragestateslocal, nonstoragestatesmean, detailedrescopl, enekvglobaldict, varendperiod)

    # Results
    updateareaprices!(price, clearing, clearing.horizons[powerhorizonix], tnormal)
    ix = vcat(ix, collect(price["steprange"]))
    pricematrix = vcat(pricematrix, price["matrix"])
    statematrix = hcat(statematrix, collect(values(startstates)))

    prices, rhstermvalues, production, consumption, hydrolevels, batterylevels = update_results(clearing, prices, rhstermvalues, production, consumption, hydrolevels, batterylevels, powerbalances, rhsterms, plants, plantbalances, plantarrows, demands, demandbalances, demandarrows, hydrostorages, batterystorages, clearingobjects, cnpp, cnph, cpdp, tnormal);    

    # Increment step
    step += 1
end

display(string("The simulation took: ", totaltime/60, " minutes"))
display(string("Time usage per timestep: ", totaltime/steps, " seconds"))

FixedDataTwoTime(DateTime("2024-12-30T00:00:00"), DateTime("1980-12-31T00:00:00"))

  5.064821 seconds (183.66 k allocations: 10.678 MiB, 1.05% compilation time)
  1.899553 seconds (2.27 M allocations: 70.734 MiB, 4.35% compilation time)
 16.189062 seconds (5.11 M allocations: 360.168 MiB, 1.54% compilation time)

FixedDataTwoTime(DateTime("2024-12-30T00:00:00"), DateTime("1981-01-02T00:00:00"))

  6.098245 seconds (72.94 k allocations: 5.026 MiB)
  0.636664 seconds (2.02 M allocations: 57.079 MiB)
 18.583206 seconds (4.94 M allocations: 351.340 MiB)

FixedDataTwoTime(DateTime("2024-12-30T00:00:00"), DateTime("1981-01-04T00:00:00"))

  4.668130 seconds (72.90 k allocations: 5.025 MiB)
  0.782235 seconds (2.02 M allocations: 57.076 MiB)
 15.299216 seconds (4.94 M allocations: 351.340 MiB)

FixedDataTwoTime(DateTime("2024-12-30T00:00:00"), DateTime("1981-01-06T00:00:00"))

  7.346379 seconds (72.92 k allocations: 5.022 MiB)
  2.614726 seconds (2.02 M allocations: 57.071 MiB)
 17.687297 seconds (4.94 M allocations: 351.412 MiB, 1.26% gc time)

FixedDataTwoTime(DateTime("2024-12-30T00:00:00"), DateTime("1981-01-08T00:00:00"))

  4.175312 seconds (72.90 k allocations: 5.022 MiB)
  0.634420 seconds (2.02 M allocations: 57.079 MiB)
  9.802127 seconds (4.94 M allocations: 351.412 MiB, 2.29% gc time)


# Only keep rhsterms that have at least one value (TODO: Do the same for sypply and demands)
rhstermtotals = dropdims(sum(rhstermvalues,dims=1),dims=1)
rhstermsupplyidx = []
rhstermdemandidx = []

for k in 1:length(rhsterms)
    if rhstermtotals[k] > 0
        push!(rhstermsupplyidx, k)
    elseif rhstermtotals[k] < 0
        push!(rhstermdemandidx, k)
    end
end

# Put rhsterms together with supplies and demands
rhstermsupplyvalues = rhstermvalues[:,rhstermsupplyidx]
rhstermdemandvalues = rhstermvalues[:,rhstermdemandidx]*-1

rhstermsupplynames = [getinstancename(rhsterm) for rhsterm in rhsterms[rhstermsupplyidx]]
rhstermsupplybalancenames = [split(getinstancename(r), "PowerBalance_")[2] for r in rhstermbalances[rhstermsupplyidx]]
rhstermdemandnames = [getinstancename(rhsterm) for rhsterm in rhsterms[rhstermdemandidx]]
rhstermdemandbalancenames = [split(getinstancename(r), "PowerBalance_")[2] for r in rhstermbalances[rhstermdemandidx]]

supplynames = [[getinstancename(plant) for plant in plants];rhstermsupplynames]
supplybalancenames = [[split(getinstancename(p), "PowerBalance_")[2] for p in plantbalances];rhstermsupplybalancenames]
supplyvalues = hcat(production,rhstermsupplyvalues)

demandnames = [[getinstancename(demand) for demand in demands];rhstermdemandnames]
demandbalancenames = [[split(getinstancename(p), "PowerBalance_")[2] for p in demandbalances];rhstermdemandbalancenames]
demandvalues = hcat(consumption, rhstermdemandvalues)

# Prepare for plotting results
hydronames = [getinstancename(hydro) for hydro in hydrostorages]
batterynames = [getinstancename(battery) for battery in batterystorages]
powerbalancenames = [split(getinstancename(getid(powerbalance)), "PowerBalance_")[2] for powerbalance in powerbalances]

# Convert reservoir filling to TWh
hydrolevels1 = copy(hydrolevels)
for (i,hydroname) in enumerate(hydronames)
    if haskey(getbalance(clearingobjects[hydrostorages[i]]).metadata, GLOBALENEQKEY)
        hydrolevels1[:,i] = hydrolevels1[:,i]*getbalance(clearingobjects[hydrostorages[i]]).metadata[GLOBALENEQKEY]
    end
end

# Time
x1 = [getisoyearstart(scenarioyearstart) + cpdp*(t-1) for t in 1:first(size(supplyvalues))] # power/load resolution
x2 = [getisoyearstart(scenarioyearstart) + cpdh*(t-1) for t in 1:first(size(hydrolevels))]; # reservoir resolution


# # Store results with binary h5 format
# datetimeformat = "yyyy-mm-ddTHH:MM:SS"
# modelname = JulES_2030_81_10

# data = Dict()
# data["areanames"] = areanames |> Vector{String}
# data["pricematrix"] = pricematrix
# data["priceindex"] = Dates.format.(ix, datetimeformat) # not necessary to store as string

# data["statenames"] = statenames
# data["statematrix"] = permutedims(statematrix)
# data["stateindex"] =  Dates.format.([getisoyearstart(scenarioyearstart) + Day(2*(i-1)) for i in 1:first(size(transpose(statematrix)))], datetimeformat)

# data["supplyvalues"] = supplyvalues
# data["supplynames"] = supplynames
# data["supplybalancenames"] = supplybalancenames

# using FileIO, HDF5
# @time h5open("$modelname.h5", "w") do file
#     for (k,v) in data
#         display(k)
#         write(file, k, v)
#     end
# end

# # # Read results so that they can be compared to other models (we have another script that compares JulES with Samnett and TheMA)
# # # Read using JLD2
# # using JLD2
# # @time JulESdata = JLD2.load("$modelname.h5")

# # # Read using HDF5
# # JulESdata = Dict{String, Any}()
# # @time c = h5open("$modelname.h5", "r") do file
# #     for key in keys(data)
# #         JulESdata[key] = read(file, key)
# #     end
# # end

# # Store as CSV
# areaprices = rename!(DataFrame(prices, :auto),powerbalancenames)
# areaprices[!,:time] = x1
# CSV.write("JulES_results\\prices.csv", areaprices)

# demand = rename!(DataFrame(demandvalues, :auto),demandnames)
# demand[!,:time] = x1
# demand = stack(demand,Not(:time))
# demandcopl = DataFrame(variable=demandnames, area=demandbalancenames)
# demand = leftjoin(demand, demandcopl, on=:variable)
# CSV.write("JulES_results\\demand.csv", demand)

# supply = rename!(DataFrame(supplyvalues, :auto),supplynames)
# supply[!,:time] = x1
# supply = stack(supply,Not(:time))
# supplycopl = DataFrame(variable=supplynames, area=supplybalancenames)
# supply = leftjoin(supply, supplycopl, on=:variable)
# CSV.write("JulES_results\\supply.csv", supply)

# hydro = rename!(DataFrame(hydrolevels, :auto),hydronames)
# hydro[!,:time] = x2
# CSV.write("JulES_results\\hydro.csv", hydro);


# Plot prices
idxwohub = findall(x -> !occursin("HUB", x), powerbalancenames) # remove hubs, not active in 2025 dataset
display(plot(x1, prices[:,idxwohub], labels=reshape(powerbalancenames[idxwohub],1,length(powerbalancenames[idxwohub])), size=(800,500), title="Prices", ylabel="€/MWh"))

# # Plot supplies and demands
# supplychart = plot(x1, supplyvalues,labels=reshape(supplynames,1,length(supplynames)),title="Supply", ylabel = "GWh/h")
# demandchart = plot(x1, demandvalues,labels=reshape(demandnames,1,length(demandnames)),title="Demand", ylabel = "GWh/h")
supplychart = areaplot(x1, sum(supplyvalues,dims=2),title="Supply", ylabel = "GWh/h")
demandchart = areaplot(x1, sum(demandvalues,dims=2),title="Demand", ylabel = "GWh/h")
display(plot([supplychart,demandchart]...,layout=(1,2),size=(800,500)))

# Plot storages
# display(areaplot(x2, hydrolevels1,labels=reshape(hydronames,1,length(hydronames)),size=(800,500),title="Reservoir levels", ylabel = "TWh")) #
display(areaplot(x2, dropdims(sum(hydrolevels,dims=2),dims=2),labels="Total reservoirs",size=(800,500),title="Reservoir levels", ylabel = "TWh")) #

display(areaplot(x1, dropdims(sum(batterylevels,dims=2),dims=2),labels="Total batteries",size=(800,500),title="Battery levels", ylabel = "GWh")) #

# Plot list of yearly mean production and demand for each supply/demand
meandemand = dropdims(mean(demandvalues,dims=1),dims=1)
meanproduction = dropdims(mean(supplyvalues,dims=1),dims=1)
supplydf = sort(DataFrame(Supplyname = supplynames, Yearly_supply_TWh = meanproduction*8.76),[:Yearly_supply_TWh], rev = true)
demanddf = sort(DataFrame(Demandname = demandnames, Yearly_demand_TWh = meandemand*8.76),[:Yearly_demand_TWh], rev = true)
supplydf[!,:ID] = collect(1:length(supplynames))
demanddf[!,:ID] = collect(1:length(demandnames))
joineddf = select!(outerjoin(supplydf,demanddf;on=:ID),Not(:ID))
# show(joineddf,allcols=true, allrows=true, nosubheader = true)

# Check that total supply equals total demand
show(combine(joineddf, [:Yearly_supply_TWh, :Yearly_demand_TWh] .=> sum∘skipmissing), nosubheader = true)

# # Plot list of yearly income and cost for each supply/demand (only works if exogenprices are collected)
# supplyrev = copy(supplyvalues)
# for (i,supplybalancename) in enumerate(supplybalancenames)
#     idx = findfirst(isequal(supplybalancename), powerbalancenames)
#     supplyrev[:,i] .= supplyrev[:,i] .* prices[:,idx]
# end
# demandrev = copy(demandvalues)
# for (i,demandbalancename) in enumerate(demandbalancenames)
#     idx = findfirst(isequal(demandbalancename), powerbalancenames)
#     demandrev[:,i] .= demandrev[:,i] .* prices[:,idx]
# end
# meandemandrev = dropdims(mean(demandrev,dims=1),dims=1)
# meanproductionrev = dropdims(mean(supplyrev,dims=1),dims=1)
# supplyrevdf = sort(DataFrame(Supplyname = supplynames, Yearly_rev_mill€ = meanproductionrev*8.76),[:Yearly_rev_mill€], rev = true)
# demandrevdf = sort(DataFrame(Demandname = demandnames, Yearly_cost_mill€ = meandemandrev*8.76),[:Yearly_cost_mill€], rev = true)
# supplyrevdf[!,:ID] = collect(1:length(supplynames))
# demandrevdf[!,:ID] = collect(1:length(demandnames))
# joinedrevdf = select!(outerjoin(supplyrevdf,demandrevdf;on=:ID),Not(:ID))
# # show(joinedrevdf,allcols=true, allrows=true, nosubheader = true)

# # Sum revenues and cost
# show(combine(joinedrevdf, [:Yearly_rev_mill€, :Yearly_cost_mill€] .=> sum∘skipmissing), nosubheader = true)

1×2 DataFrame
 Row │ Yearly_supply_TWh_sum_skipmissing  Yearly_demand_TWh_sum_skipmissing 
─────┼──────────────────────────────────────────────────────────────────────
   1 │                           4226.69                            4221.37

Prototype JulES v1¶

Import packages¶

Prepare parallell processing - import code on all cores¶

Starting point and scenario settings¶

Price prognosis models uses NVEs Europe dataset for the TheMA model¶

NB!¶

Initialize and run price prognosis models for all scenarios¶

Plot medium term prices for some scenarios¶

Plot short term prices for some scenarios¶

Stochastic subsystem models and market clearing model uses more detailed hydropower data¶

Get water values for detailed hydro (DETD) from aggregated hydro (TheMA)¶

Scenario modelling¶

Stochastic subsystem models with Benders decomposition¶

Market clearing problem¶

Head dependency in clearing_init!() and clearing!()¶

Initialize and collect results: Prices and start states (e.g. reservoirs)¶

Collect more detailed results¶

Simulate next time steps (16 steps = 32 days)¶

Results¶

Postprocess detailed results¶